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Abstract 

In a series of publications the Cummings-Stell model (CSM), for a binary mixture of asso- 
ciative fluids with steric effects, has been solved analytically using the Percus-Yevick approx- 
imation (PYA). The solution consists in a square well potential of width w, whose center is 
placed into the hard sphere shell (r < a): at L — a/n (n — 1, ... ,4). This paper presents 
a general solution, for any n, of the first order Difference Differential Equation (DDE), for 
the auxiliary Baxter's function that appears in the CSM, using recursive properties of these 
auxiliary functions and a matrix composed by differential and shift operators (MDSO). This 
problem is common in some other models of associative fluids such as the CSM for homogeneus 
and inhomogeneus mixtures of sticky shielded hard spheres including solvent effects under PYA, 
and in that of mean-spherical approximation (MSA), for chemical ion association and dipolar 
dumbbells and polymers. The sticky potential implies a discontinuity step at L in the solution 
of auxiliary Baxter's functions so that, one side, L now is arbitrary and, for some additional 
effects, it can be placed one or more sticky potentials at different positions into the hard shell. 
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1 Introduction 

In 1984 Cummings and Stell [Ij proposed a simplified haniiltonian model for a reactive system of 
two types of homogeneus fluids with the same density and diametei0- In this model, type A and 
B molecules, can be associated by means of a selective square interaction that, in the appropiate 
limits, can be reduced to a sticky Baxter's potential [2 located inside the hard core at a distance 
from the center of the particle L = a/n. The cases of n = 2, 3, 4 for the Baxter's function have been 
solved in [ll[3l[4] using the same formalism and, at [5], they apply the results for a pair of reactive 
fluids and a solvent. Lee and Rasaiah solved the L = a/A, L = a/5 and proposed a solution for 
a/n for chemical association and dipolar dumbbells [6]. 

The model of Lee and Rasaiah of association in electrolytes + ^ AB is studied in [7] . In 
this case the authors add a selective coulombian part, ±e^/r, out of core, to the associating original 
CSM. 

The sticky site inside the hard core incorporates geometrical conditions of steric saturation in 
the molecule and this idea is shown, using computational simulations in different ensembles, in 
ref. [8]. For different bonding length parameters the system allows formation of dimers for small 
L, chains for L slightly larger, and vulcanization of species for bonding length values close to the 
diameter a of particles. Huerta and Naumis studied the connectivity of a binary mixture using a 
selectively hard sphere potential and a superposition as [9J: 

U.jir)^U^^{r) + il~S.,)UUr) 

where 

" ^ ' 10, r > 1 

r < L — 0.5w 
L - 0.5w < r <1 
r > 1 

r < L — 0.5w 
Uas{r) = { -Eas ~ D, L- 0.5w <r <l + 0.5w 

r > L + 0.5w, 

where L is the bonding distance, w the intracore square well width and i, j represents the species 
in the mixture. The flnal potential (see Figure [T]) is equivalent to the original Cummings-Stell in 
the adecuate limits for sticky approximation. 

The same idea is implemented by Pizio and Blum |10j for a hard-sphere fluid with dimerization 
A + A ^ A2. In the development most of the models maintain iy as a parameter (bonding distance) 
and flnally take the case L = a/2, however some other possibilities are presents having the analytical 
solution for arbitrary L. Kalyuzhnyi and Stell [llj present a recount of cases for different ranges 
of the location L. As we show here, the sticky potential into the hard sphere shell produces a 
discontinuity step in the auxiliary functions of Baxter. This fact allows us to think of systems with 
more than one sticky site inside the shell or, even, a distribution of sticky wells. 



and 




^It will be detailed in next section. 
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Figure 1: The image shows the basic form of the potential, 
different L according to [9] . 



Below, their structural effects for 



The vertical dashed arrows show the sense of the Baxter's sticky limit pQ. The molecular 
diameter is a and the dashed curves (with labels "AA" and "AB" ) corresponds to the colulombian 
interaction for electrolyte [?]■ 

In the following sections, we first develop the matrix of differential and shift operators MDSO, 
and its inverse, for the simplest case n = 2, followed for the inversion of the MDSO's for a very 
general case that corresponds to the set of n DDE's, to n subintervals of [0, a] and to a sticky 
location L = ma/n. 
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2 The association model: hard spheres with shielded sticky 
interaction 



2.1 The model of binary mixture 

The statistical mechanical model of chemical reactions of Cummings and Stell [T] represents the 
association of two species A + B ^ AB, with the same density and diameter, which simplifies the 
mathematical problem. The potential proposed in the CSM consists of a hard sphere repulsion 
between like species (A-A or B-B) and a mound of width a with a deep, narrow, and attractive 
square well, with width w centred on L. Here L < cr/2 and L + w/2 ^ a/2 for AB interactions: 

ei ifO<r<L-w/2 

' ' ei \iL + w/2<r<a 

iir > a 

The geometric consideration of this model for the AB interactions ensures steric saturation in the 
system (there is no formation of 7i-mers for n ^ 3) due to overlapping. In addition, this model has 
a solution in the PY approximation, mapping the square well onto an infinitely deep and stretch 
well like the sticky potential of Baxter [2 . The connection between the Baxter's original and this 
model, is obtained equating the second virial coefficients: first it is considered the limit ei — > oo, 
€2 — > oo, w ^ 0. The limits of £2 — > oo and w — > are taken to maintain tractable the problem in 
the PY approximation [2]. The limit ei — > 00 (in the repulsive part: see figure [1} doesn't change 
very much the results, but simplifies the solution [3]. 

The total and direct correlation functions are related by the Ornstein-Zernike (OZ) equation 
that, for this binary mixture, can be written as [12j 



/iy(r) = Cy(r) + Pk,o Cik{s)h^^{\r - s\)ds (2) 



k=A,B 



where the integral is evaluated in the whole space, s — \\s\\, and pkfi is the number density of 
species k particles. In the CSM the densities are considered equal. 

Briefly reviewing the formulation of the CS model: we need the factorized form of the OZ ^13j 
equations ([2]) with the PY closures 



hAA{r) ^ hBB{r) = -I, r < a 
CAA{r) = CBB{r) = 0, r> (7 



(3) 



and, in the considered limits 

hAB{r) = -l + ^5{r-L), r < a \ 

(4) 



CAB{r) = 0, r> a. 

Given the conditions of the problem, the factorized OZ equations are written [Ml [U [3l |4] 

rhii{r) = -<zii(r) + 2^p dt{r - t) [gn(t)/tn(|r - t\) + <Zi2(t)/^i2(k - 

(5) 

rhi2{r) - -g'laW + /o" d^ir - t) [gn(t)/ii2(|r - t\) + qi2(t)/iii(|r - t\)] . 
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where we changed the index of Haa or Hab to hn or hi2 respectively. The same for ah functions. 

Substituting the closure relations (|3]) and (0]) in the set of OZ equations ([5]), the following system 
of difference-differential equations (DDE) for the auxiliary Baxter's functions qij (r) is obtained: 



9n ('") + p[Qi2 {r + L) - qi2 {r - L)] = (an + Dai2)r + bu + Dbi2 

(6) 

q[2ir) + P[qii{r + L) - qii{r - L)] = {Dan + ai2)r + Dbn + bu - ^5{r - L), 



where p — irpXL /6 and 



Sij - 211 p qij{t)dt 



hj = 2'Kp tqij{t)dt 
and satisfies the boundary conditions 

qii(cr) = (7i2(o-) = 
qi2{L-) = qi2{L+) + ^ 



(7) 



After integration, a step appears in r = L for the auxiliary function qi2{'r) due to the delta term 
associated with the well. 



2.2 Solving for n = 2 

In the first work of Cummings and Stell [T] a new pair of functions were used and defined as the 
sum and difference of the originals qii(r) and (712 (r). The advantage of this trick is to obtain two 
uncoupled equations, one for q+{r) and another for q-{r) which can be solved in a separate way. 
If the functions q+{r) — qii{r) + qi2{'i') and q-{r) = qii{r) — qi2{r) are defined then, adding and 
substracting Eqs. ([6]), 



q+{r) +p[q+{r + L) - q+{r - L)] = a+r + b^ 
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S{r - L) 



and 



q'_ (r) - p[q- 
with the obvious definitions^ 



and 



+ L) - q^{r - L)] = a_ 



(8) 



(9) 



2ttp I dtq±{t) 
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6± = (1 ± I?)27rp / dttq±{t). 
Jo 



The last term in ^ and Q can be omitted since it is equal to zero in all subintervals except where 
r = L ~ ajn. This condition is fixed in the boundary conditions ([7]). In the rest of subintervals, 
for the general case, must be true that [1] [6] 



qimajn ) = (/(mci/n^), for m — 2, 3, . 



1. 



(10) 



^From [4] D = 1, so that a_ 
paper. 



: 0. This not imphes changes in the results. We asume this fact in the rest of 
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We use here q{r) = qii{r) + qi2{r). The aim of this proposal is to find an analytical form for 
the function q(r) assuming that: i) the solution must be made in subintervals |15| . ii) this implies 
that g(r) will be defined also in subintervals, and iii) the original functions qij{r) can be recovered: 
= ^'^^''^2*^'^''^ and qi2{r) = . Identical procedure shows that is sufficient to replace 

A —>■ —X,p — » — p, ly —>■ —V to obtain q~{r). 



3 The MDSO 

The cases shown in [TJ [3J HJ |S] are solved here using MDSO. For convenience we show the case 
L = (t/2 of the CS model [T] in detail, and the cases L — u and L — a/A summarized. 



3.1 The case n = 2 

The first case yields the system of coupled differential equations 

dqi{r) 



dr 

and 



+ pq2{r + cr/2) =ar + &, forO < r < a/2 (11) 



^^^^ - Ml (r- cr/2) = ar + fe, for ct/2 < r < cr. (12) 
dr 

where, evidently, qi corresponds to the first half of the interval and q2 to the second. We define 
here the differential operator V as 'Df{x) = ^^^^ and the shift operator 8^ by f{x) = f{x ± s). 
With this operators defined, the set of (fTTj) and p2)l can be rewritten as 



and 

or, in matricial form, as 



Vqi{r) +p£''/^q2{r) =ar + b 
Vq2{r) -pS^^/^qiir) =ar + b 



V \( qi{r)\_( h{r) 

V \q2{r) -\ f2{r) 



(13) 



These equations can be reduced to a symbolic form 

M2q{r) = f{r) (14) 

where is the matrix of differential and shift operators, or MDSO, that appears in (|13p. applied 
to the vector q of functions qi{r). The right side is the vector / of functions fi{r) that, in this case, 
are linear functions of r. The index in M, corresponds to the number of equations (or partitions 
in the interval of solution). 

The main idea of this paper is to find a solution for the system represented in ()14p as 



qir)=M^'f{r). 

This implies the knowledge of an explicit analytical form of the inverse of 2 , and how it operates 
on f{r). One way of defining the inverse of the differential operator V is by using the equation 

y'{x)±ay{x)^ fix) (15) 
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{V ± a)y{x) = fix) 
whose solution leads us to define the inverse operator {V ± a)~^ as 

(V ± a)~^f{x) = Ce^"'' + e^F-^^ j e^"'"' fix')dx' . (16) 

In the previous expression, the case a = imphes that the inverse MDSO is reduced to the trivial 
definition of inverse differential operator as an integral operator. The case where a is a complex 
number (or a pure imaginary one) implies harmonic solutions [16] and Fourier transform of the 
right hand side of differential equation. 
Continuing with the case L = a/2, the inverse of A^2 is 



^2 



where the commutation properties of the operators V and £ were used. Direct calculation gives the 
determinant-operator of A^^^ a|l 

A^ = V^+p^ = {V + ip){V-ip), (18) 

so that (fT7|) and ([T5)l define completely the inverse determinant-operator of as the product of 
two inverse operators of the form of ([16]) : 



A2 V + ipV — ip 

This is the formal inverse determinant of the MDSO, however we still need to find the appropiate 
coefficients to satisfy the boundary conditions. So that, the direct application of the inverse MDSO, 
dl), on (dl]) and ^ gives [HIIIH] 

qi{r) = Acospr + Bsmpr-^r+p{l-iy/2)-^ 
q2{r) = Ccoapr + Dsmpr + ^r+^{l~iy/2) + ^ 
with p = pa. Now, considering ([7]) and the fact of p2| must be satisfied we obtairQ, explicitly. 



gi(r) = Acospr + Bsmpr- !ir+^il-iy/2)~ ^ 

q2{r) = Aainpir - a/2) + B coapir - a/2) + ^r+ ^{1-1^/2) + I 



(20) 



which agree exactly with the results in IJ . The second equation, now has the same set of constants 
that the first. The harmonic functions have been interchanged and their arguments are shifted by 
-a/2. 

•^The application of inverse shifting operators is the identity: £" S~" f{x) = f{x + s) = f{x + s — s) = f{x). 
^This fact allows us to establish the same set of constants for the harmonic part of the solution. 
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3.2 The more general case n > 2 

In the case L = a/'i there are three equation^: 



dr 

dq2{r) 



^^^'■'^^ +pq2ir + a/3) ^ ar + bJorO < r < a/3 



dr 

with the MDSO given as 



+ pq3ir + CT/i)-pqi{r - cr/3) = ar + 6, for cr/3 < r < 2cr/3 (21) 
Pq2{i' ~ cr/3) = ar + 6, for 2cr/3 < r < a 



dr 

dqsir) 



V pE"/^ 

= I -pE-"'^ V pE"/^ 
-pE-"/^ V 



and the inverse operator of Ai 
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^ / V^+p^ -pVE''/^ p^E^''/^ 

■^3^ = — pVE-""/^ P2 -pVE"^'^ 

V{V +2p ) y ^2g-2a/3 p2^p2 



where 



A3 V(V^ + 2p^) {V - 0) (V + iV2p) [V ~ tV2p) 



/ I?2+p2 _p2?£-/3 p2^2./3 \ 

= pVE-"l^ ^pVE"/^ . (22) 

V p2£-2./3 pP£-^/3 2?2 ^ p2 y 

Obviously I/A3 is a product of inverse operators in the form of Applying these to the right 

hand side of ([2T|) we obtain, directly 

qi{r) = Ai cos V2pr + Bi sinv^pr + ^r^ - ■;^(1 - 2iy/3)r + 6r + Fi 

2 2p 



<3'2 



(r) = A2 cos v^pr + B2 sin V2pr + F2 



qsir) = A3 cos V2pr + B3 sin %/2pr + ^r'^ + ^{1- 2v/3)r + br + F3 

I Zp 

and, imposing bound conditions in the respective subintervals and the fact of (j2ip must be satisfied 
(as in (PO)) for ri = 2), we obtain v42, A3, i?2 and B^ in terms of Ai and i?i 

gi (r) = Ai cos V2pr + Bi sin V2pr + ^^^^ - (1 - 2iy/3)r + br 

2 2p 



^The first and last equations always have one term less, due to the condition of PYA, q{r) = out of [0, cr], 
^With this we are defining Mn^ = -^A4n- 
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92 



(r) = V2A1 sin V2p{r - a/'S) + V2Bi cos V2p{r - ct/3) 



(23) 



qsir) = -Ai cos V2p(r - 2cr/3) - Si sin\/2p(r - 2ct/3) + ^r^ + -^(1 - 2v/3)r + br. 

2 2p 



The case i = a/ A has a tridiagonal matrix Ai^, whose determinant A4 = + SV^p^ has 
roots ±i^^^p and ±i^^^p so that 



A4 = 



V + I — - — p 



^ V5- 1 
u — I — - — p 



P + ^ p 



— « — - — p 



Lee and Rasaiah, in p], caUed these roots x = and y = -^^^3^, and correspond to the a^'s 

defined below in this paper. 

The case n = 5 or L = cr/5 has a determinant A5 — T)^ + AV^p^ + 3'Dp'^ whose roots are 
0, ±ip, ±i-\/3p, and the solutions have the same structure as ([201) and (|23p . and that showed in 



4 The general MDSO 
4.1 The matrix 

Making the same construction for n divisions in the solution interval, one obtains n functions qi 
that represent a continuous solution. Each of them must be such that 

■ 1 i 
—a < r < —a 

n 



Qii"^) +P['li+i{r + <^/n) - qi-i{r - cr/n)] ^ ar + b, for 

n 

where i = 1, 2, 3 . . . , 71 and go = 9^+1 = 0, in which 

q{r), 



0, 



-a < r < -a 

(, n 

otherwise . 



The generalized MDSO for arbitrary n has the form tridiagonal 



V 

pE"! 



for I = j 
for i = i — \ 
pg-<y/n fori = j + l 

oherwise 



or, in matrix form 



\ 



V 


p£' 








pS-" 


V 


pS" 








-pS-" 


V 


p£^ 










-pE-" 















p8^ 
V 

^p£- 



\ 






p£^ 
V ) 



(24) 



(25) 



(26) 



Due to symmetry of A^„ the shift operators are mutually canceled in the inverse of A„. This fact 
enables us to put all solutions qi{r) in terms of inverse differential operators of the form (|16p and 



9 



4.2 The determinant inverse operator 

From the tridiagonal matrix obtained, (|26p . the determinants for different values of n can be 
evaluated: 

Ao = 1 

Ai = V 

A2 = V^+p^ 

A3 = + 2Vp^ 

A4 P"* + 3V^p^ + / 

A5 = P'"^ + 4p3p2 _^ 3pp4 

Ae = I?^ + 5I?V+6PV+/ 

Ai4 = v^^ + isp^V + eep^V + i65x>V + 2101?^^ + i26v^p^° + 

+28V^p^^ +p^'^ 

Ai5 = V^^ + UV^^p^ + 781?"/ + 220X>V + SSOPV + 252V^p^" + 



It is easy to prove that the recurrence relation between determinants of MDSO's of different order 
is, for n ^ 2, 

A„==I?A„_i+/A„_2 (27) 

where Ai = T> and Aq = 1. The index of A„ corresponds to the order of the MDSO. From this 
recurrence relation we obtain the general expression for A„ 

A„ = v-+j:T=i^^ut7i^-k) 

Em pt-ajp^j („_j)! 
j=0 j\ (n-2j)! 

which, finally, can be reduced to 

n~ J 
where, additionally 

n 

An^Y[{V-Xk). (29) 

fe=i 

In these expressions m =[n/2j is the integral part of n/2. In we have written A„ in factorized 
polynomial form. Here Xk are all the n roots of A„ which are all pure imaginary and proportional 
to p. When n is even there are exactly m = n/2 pairs of complex conjugated root^ ±iakP and, if 
n is odd, there is a further null root of A„, which requires an additional integration to obtain the 
qi functions. 



A„=^( jp-V^- (28) 



^This was proved for n = 1, 2, . . . , 26. 
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We can write A„, with m as defined above, as 

nr=i^fc for n even 



(30) 



wliere we define = {T> ~ Xk){T> Xk) and Xk = isk = icv-kP- Also, Xk is the complex conjugate 
of Xk, OLk is a real number, and I?o is the operator associated with the root xo = Q for n odd. 
Direct application of individual inverse Vj. on a linear function gives [16] 



■Dj, ^(cr + d) = Ak cosskr + Bk sinsfer + — (cr + d). 



(31) 



Same as before, Sk is defined by mean of Xk = «Sfe = iakP and D^. ^ = {Df.) ^ ■ In the general case, 
a partition n even of the interval [0, a] results in 



V^^V^ ^ ... V;^{cr + d) = Y^ Ak cos s^r + Bk sin Sfcr 



fc=i 



n s2 



(cr + d). 



(32) 



It is easy to see that if we define the vectorial funcion v{r), using the definition in (fT7|) or (p2)l . as 
the application of 7Vi„ to the vector /(r), 

v(f) = Mnf{r), 

we obtain a linear function t>(rjf|. The general solution of the system of equations become 



q(r) = -^v{r) 



where 



q{r) = 



( qiir) \ 

92 (r) 
V Qn{r) J 



and qi{r) is as defined in The constants and i?i that appears in ([5^ can be stablished 

from bound conditions at the frontiers of the subintervals given by ([7]) and (jlOp and the recurrence 
stablished in the original set of equations. 



5 Discussion 

A general solution of the set of DDE's for the original CS model has been discussed here. It requires 
the roots of the polinomial expresion for the inverse determinant A^^ which always are of the form 
iskP with Sfc a real number. One of the advantages of this method is the fact that it allows to 
choose the site of the sticky potential, not only at L = a/n for a few values of n: the solution is 
valid for any n. It allows to change the sticky potential site to positions more and more close to (or 

*The harmonic part comes from the inverse determinant. 
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away of) the center of the particle. For example in the case n — 3 the step of discontinuity XL^/12 
can be defined at i = a/3 or at L = 2a/3. 

This fact bring us the possibility to set the well at r = ma/n to obtain different molecular 
structures [9]: by locating the step in the site of sticky and represent a step discontinuity 

in the solution at the position L = ma/n. With this in mind one can think, also, the posibility of 
two or more sticky square wells into the hard shell. 

Acknowledgements 

J. F. Rojas, is very grateful to Prof. L. Blum. 

References 

[I] P.T. Cummings and G. Stell. Statistical mechanical models of chemical reactions. Analitical 
solution of models A + B ^ AB in the Percus-Yevick approximation. Mol. Phys., 51:253, 
1984. 

[2] R.J. Baxter. Percus-Yevick Equation for Hard Spheres with Surface Adhesion. J.Chem.Phys., 
49:2770, 1968. 

[3] P.T. Cummings and G. Stell. Statistical mechanical models of chemical reactions II. Analitic 
solution of the Percus-Yevick approximation for a model of homogeneus association. Mol. 
Phys., 55(l):33-48, 1985. 

[4] P.T. Cummings and G. Stell. Statistical mechanical models of chemical reactions III. Solvent 
effects. Mol. Phys., 60:1315-1342, 1987. 

[5] S. H. Lee, P.T. Cummings, and G. Stell. Statistical mechanical models of chemical reactions 
IV: Solvent effects near the critical point. Mol. Phys., 62:65-90, 1987. 

[6] Song Hi Lee and Jayendran C. Rasaiah. Chemical ion association and dipolar dumbbells in 
the mean spherical approximation. J. Chem. Phys., 86(2):983-994, 1987. 

[7] Song Hi Lee and Jayendran C. Rasaiah. A model for association in electrolytes, analytic 
solution of the hypernetted-chain/mean spherical approximation. J. Chem. Phys., 83(1):317- 
325, July 1985. 

[8] A. Huerta. A two-dimensional model associating fluid with spherically integral equations and 
monte carlo simulation study symmetric intracore square well shell. Mol. Phys., 96(5):795-804, 
March 1999. 

[9] A. Huerta and G.G. Naumis. Relationship between glass transition and rigidity in a binary 
associative fluid. Phys. Lett. A, (299):660-665, 2002. 

[10] Orest Pizio and Lesser Blum. Analitic solution of the mean spherical approximation for a dipo- 
lar hard-sphere fluid with intracore anisotropic sticky interactions. Phys. Rev. E, 52(l):572-579, 
July 1995. 

[II] Yu V. Kalyuzhnyi and George Stell. On the effects of association in fluids with spherically 
symmetric interactions i. Mol. Phys., 78(5):1247-1258, 1993. 



12 



[12] L. S. Ornstein and F Zernike. Accidental deviations of density and opalescence at the critical 
point of a single substance. Proc. K. Ned. Akad. Wet, 17:793-806, 1914. 

[13] L. Blum, J. F. Rojas, and J. N. Herrera. Correlation functions I: formalism and simple appli- 
cations. Rev. Mex. Fis., 39(5):799-817, 1993. 

[14] Peter T. Cummings. Analytic solution of integral equations for molecular fluids. Kinam, 
16(Scric A):121 141, 1984. 

[15] R. Bellman and K.L. Cooke. Differential- Difference Equations. Academic Press, 1963. 

[16] Guadalupe M. Munguia Gamez and Martin G. Garcia Alvarado. Ecuaciones diferenciales con 
retardo. Memorias de la XVII Semana Regional de Nivel Superior Investigacion y Docencia 
en Matemdticas. Mosaicos Matemdticos., (20), August 2007. 



13 



